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Abstract. - Employing a recently developed dynamical density functional theory we study the 
response of a colloidal sediment above a wall to shear, demonstrating the time dependent changes 
of the density distribution and its center-of-mass after switching shear either on or off and under 
oscillatory shear. Following the onset of steady shear we identify two dynamical mechanisms, 
distinguished by their timescales. Shortly after the onset, a transient enhancement of the packing 
structure at the wall reflects the self-organization into lanes. On a much longer timescale these 
effects are transmitted to the bulk, leading to migration away from the wall and an increase 
in the center-of-mass. Under oscillatory shear flow the center-of-mass enters a stationary state, 
reminiscent of a driven damped oscillator. 
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Colloidal suspensions under shear have been the sub- 
ject of intense research during the last decades, both the- 
oretically as well as in simulations and experiments [TfH]. 
Dense systems have attracted particular attention, as par- 
ticle interactions can transmit the effects of shear into di- 
rections perpendicular to the flow, leading to nontrivial 
collective dynamics. An important consequence of apply- 
ing steady shear is that the intrinsic slow dynamics ex- 
hibited by quiescent dense systems (and characterized by 
observables such as density correlators or mean squared 
displacements [5H9]) can be sped up drastically in all spa- 
tial directions, provided the shear rate exceeds the in- 
' verse of the structural relaxation time [10| . The dynam- 
ical response becomes richer still when considering time- 
dependent shear fields for which transient dynamics can 
play a dominant role (e.g. oscillatory shear [XT]). 

In real colloidal systems the buoyant mass of the sus- 
pended particles is often nonzero, such that gravity drives 
the formation of a nonuniform density distribution: A col- 
loidal sediment (or creaming profile, in the event that the 
buoyant mass is negative). The study of colloidal sedi- 
mentation has a long history, effectively starting with the 
pioneering experimental work of Perrin at the start of the 
19th century. While early theoretical works focused upon 
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the influence of hydrodynamic interactions on the sedi- 
mentation velocity in dilute systems (see e.g. [T2H14] ). the 
dynamic settling of denser suspensions has only recently 
been studied in detail in experiments, Brownian dynam- 
ics simulations, and theoretically, using dynamic density 
functional theory (DDFT) [T5h17| . Incorporating hydro- 
dynamic interactions into theories of dense suspensions is 
a difficult task, although recently some progress has been 
made [HJQS]. 

Experiments in which shear flow has been applied to 
particle sediments (where in the early literature, mostly 
large, nearly non-Brownian particles were considered) re- 
port changes in the viscosity, attributed to a shear in- 
duced modification of the underlying microstructure |20j . 
Moreover, an increase in the height of the sediment has 
been observed with increasing flow rate. These two obser- 
vations imply that shear flow acts to reduce the density 
within the sediment and redisperse the particles. This ro- 
bust physical effect, now commonly referred to as 'viscous 
resuspension', has been addressed in a number of studies 
20 24! and has clear relevance for commercial and indus- 
trial systems in which unwanted sedimentation effects can 
be suppressed by mechanical means. 

In this paper, we use the framework of DDFT to study 
the interplay of sedimentation and shear, concentrating 
on situations for which the direction of shear flow and the 
gravitational force are perpendicular. Standard DDFT 
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[2"5] for flow advection does not capture the effects of shear 
in this situation. It lacks the essential physical ingredient; 
namely the transmission of shear effects into directions 
perpendicular to the flow. We therefore employ a recently 
suggested modification of the theory which corrects this 
failing [25]. By considering the time dependence of the 
particle distribution we obtain insight into the microscopic 
mechanisms at work during the process of viscous resus- 
pension. 

Following the onset of shear we observe a two step sce- 
nario. On a short timescale, particles close to the wall 
build lanes to allow for lateral movements following the 
shear flow. This local self-organization leads to a tempo- 
rary reduction in the entropy of the system and drives the 
slow migration of particles to regions far away from the 
wall. Only on a much longer timescale do we observe an 
increase of the center-of-mass (COM) and, therefore, in 
the gravitational potential energy, leading eventually to a 
steady state. If the shear field is then suddenly switched 
off, we find that the equilibration dynamics show an inter- 
esting symmetry with that following switch on, at least for 
moderate shear rates. In the case of oscillatory shear the 
transient regime occurring in response to switching on the 
flow is followed by a periodic stationary state, out of phase 
with the applied strain. For large frequency the COM con- 
verges to a time independent value above the equilibrium 
value. Given the possibility of precise real time particle 
observation and tracking in shear cells, e.g. using confocal 
microscopy [27] , we hope that our theoretical findings will 
motivate further experimental studies of the dynamics in 
sediments subject to shear flow. 

Consider a sediment of N spherical colloidal particles, 
i = 1 . . . N, of diameter d, dispersed in an incompressible 
Newtonian solvent. The time evolution of the probabil- 
ity distribution of particle positions, ^(t) = ^({r^}, t), is 
given by the Smoluchowski equation [1] 



at 



0, 



where the probability flux of particle i is given by 

j;= [/c(t)Ti-IV(0i-/3Fi)]*(t), 



(1) 



(2) 



where /? = (^T)" 1 is the inverse thermal energy and Dq 
is the (bare) one particle diffusion coefficient. n(t) is the 
velocity gradient tensor, which we choose to be j(t)xy, 
i.e., we consider time dependent simple shear pointing in 
x direction and varying in y direction. Introducing dimen- 
sionlcss units gives space, energy and time in terms of d, 
ksT and d 2 /D , respectively. The shear rate 7 is in these 
units equal to the Peclet number Pe= jd 2 /D . 

The force on a given particle Fj = ~djU is generated 
from the total potential U, made up of a pairwisc additive 
interaction <j> and an external potential V ext , 



Udn}, t) = £ V* ({r,}, t) + £ 0(|r< - r 3 -|). 



(3) 



The sediment of colloidal particles shall be confined to 
y > 0, i.e., we introduce a repulsive wall in the plane 
y — 0. Gravity acts perpendicular to the wall, i.e., in 
direction -e^, and the external potential 



V c * t ({r l },t) = V c * t ({y l }) = 



mgy l 



(4) 



depends on the y component of the particle positions only. 
The Yukawa potential (shifted by -1/2, as the center of 
the spheres is restricted to y > 1/2) carries a parameter s 
controlling the softness of the repulsion, and mgyi is the 
gravitational potential of particle i with buoyant mass m. 
The pair potential <fi is taken to be a hard sphere poten- 
tial, which well describes the interaction between certain 
colloids, 

<t>{r) = { ° r> ! • (5) 
^ w [ 00 r < 1 w 

The probability flux ([2]) neglects hydrodynamic interac- 
tions and assumes that the velocity gradient tensor n(t) is 
translationally invariant in space. We refer the reader to 
a detailed discussion of these approximations in Ref. [4] . 

We use dynamic density functional theory (DDFT) to 
calculate the one-body density profile, which arises from 
"J by integration over N — 1 particle positions [T] 



p(r,t)=N 



T% ■ 



d 3 rjv*({rj,t). 



(6) 



For the potential in Eq. ((3]), one has p(r,t) = p(y,t) [26] . 
As discussed in detail in Ref. [55], standard flow advected 
DDFT [25] does not capture the effects of shear in the 
considered setup, because the flow advection term van- 
ishes exactly: It does not transmit the effects of shear, 
pointing in x direction, to particle motion against the ex- 
ternal potential in y direction. It is therefore necessary 
to introduce a modified version of the DDFT equation for 
sheared systems [26] 



dp{y,t) = d 
dt dy 

The first term on the right hand side captures the velocity 
of particles in the y direction due to shear in a mean field 
manner. As a true microscopic derivation of the shear 
term is still lacking, it must be viewed as an ad hoc, albeit 
physically motivated, addition to the theory [2l)] . 



vl{y,t)= / dy'p(y',t)v*(y-y',t). 



(8) 



Eq. © is a convolution of the density with a velocity Vy, 
describing the collision process of two particles with sep- 
aration y — y' . Within the simplest approximation it con- 
tains the equilibrium pair correlation function g cq (d), 
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(9) 
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Fig. 1: Evolution of a colloidal sediment above a Yukawa wall 
following the onset of steady shear flow (7 = 10). Inset (a) 
shows the height of the first peak as a function of time for 
7 = 10 (considered in the main figure) 7 = 15 and 7 = 20. 
Inset (b) shows the long range decay of the profile at r = 
and 00 (full black curves), and at intermediate times r = 1, 5 
and 10 (broken red curves). 

Self-consistent updating of g eq {d) to make it shear depen- 
dent as well as incorporation of hydrodynamic interactions 
are possible, in principle. The last term in (J7]) is derived 
from the Hclmholtz free energy 

T[ P (y,t)] = T c *{p(y,t)} 

+ J d 3 rp(y, t) [ln(p(y, t)A 3 ) - 1 + V c ^(y)] , (10) 

where the excess free energy 7 7CX depends on the inter- 
action potential </>. For our hard sphere system, we use 
the well known Rosenfeld approximation [28]. A is the 
thermal de Broglie wavelength. The starting point for our 
dynamical calculations is an initial equilibrium sedimen- 
tation profile with adsorption T = 3.7864 (equal to the 
total number of particles in a vertical column of unit cross 
sectional area) calculated using the parameters s = 5 and 
mg = 0.7. The average particle number N (albeit formally 
infinite as there are no boundaries in x and z directions) 
is fixed. 

We first consider discontinuously changing the shear 
rate at t = t on from zero to a constant value 7 = 10 
[27,29], denoting r = t — t on . FigfT] shows the correspond- 
ing time evolution of the density. Immediately following 
the onset of shear, the density increases at the peaks and 
reduces between the peaks, i.e., the peaks sharpen: This 
laning in the vicinity of the wall minimizes collisions and 
allows the particles to move past each other. It is only on 
much longer times that the height of the first peak reduces 
and the tail of the profile grows monotonically to larger 
y values, saturating to the steady state solution for strain 
7T = 0(100). 

Inset (a) shows the nonmonotonic behavior of the height 
of the first peak (all other peaks evolve in a similar fashion) 



as a function of time, both for the shear rate considered 
in the main figure (7 = 10), as well as for two additional 
higher rates, 7 = 15 and 20. The time at which the max- 
imum of the first peak overshoot occurs is not a sensitive 
function of shear rate (maxima occur at r = 0.240,0.251 
and 0.272 for 7 = 10,15 and 20, respectively). For the 
cases shown, the maximum height of the overshoot in- 
creases slightly faster than linearly with 7. Inset (b) fo- 
cuses on the tail of the profile, which gradually increases 
with time, reflecting the slow collective diffusion of par- 
ticles away from the wall. During this slow process, the 
curves in inset a) approach their final values close to the 
equilibrium value. The general similarity of the first peak 
for r = and r — > 00, both in height and form, is related 
to the fact that the overall weight of particles supported 
by the wall is equal in both situations (An analogous cal- 
culation for a strict hard- wall external field yields identical 
values for the contact value p(R,r < 0) = p(R,T — > 00) 
[26]). 

A clear difference between steady and equilibrium pro- 
files is seen in the vicinity of the wall (1 < j/ < 3), where 
the former is clearly lower. This is a signature of particles 
having moved away from the wall into the tail of the dis- 
tribution by collision induced diffusion (the total particle 
number is conserved). 

The "lifting" of the density distribution against grav- 
ity after switch on can be quantified by the COM of the 
density distribution 
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FigEk shows the COM as a function of time after switch- 
on. We observe the expected monotonic increase, starting 
at strains greater than unity and saturating to a steady 
state value at strains ~ 300. In [26] we demonstrated that 
the final value of the COM increases with shear rate. The 
change of the density near the wall in Fig. Q] leads to a 
change of the osmotic pressure acting on the wall, given 
by 



n w (t) 



dy' P (y',t) dV ^} y) (12) 



(which reduces to the familiar sum-rule /3II = p(0 + ) for 
the case of a hard wall). For both r < and r — » 00, II w (i) 
equals the equilibrium pressure TL eq — mg T. In Fig. [2b 
we show the excess pressure n ex (i) = n w (i) — n cq as a 
function of time. This excess pressure lifts the particles 
against the Stokesian friction and must equal the average 
drag force. In the figure, we verify numerically that 



n, 



j_ r dy(t) 
jjL at 



(13) 



with p the mobility (included here explicitly, despite our 
dimensionless description) is fulfilled at all times. 

Noting the logarithmic timescale, the excess pressure 
also shows the range of timescales involved, reaching a 
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Fig. 2: (a) The time development of the COM following the 
onset of steady shear flow for 7 = 10. Saturation to a con- 
stant value requires strains 7T in excess of 100. (b) Excess 
osmotic pressure acting on the wall. The line is the result 
calculated from Eq. after subtraction of the equilibrium 



using Eg . (|13fl . The osmotic pressure has a maximum at a strain 
O(l). (c) The free energy gain of the system as a function of 
strain. The dashed line shows the loss of free energy after 
switch off, see main text, (d) The entropy change as a function 
of strain. The dashed line is AS* = and serves as a guide to 
the eye. For all curves we employ the same parameters as used 
in FiglU 



maximum at t — 0.213 (which roughly coincides with the 
maximum of the first peak in Fig. [TJ 7T = 2.4), before 
decaying over a much longer range of strain (jr < 1000), 
corresponding to the timescale where the particles migrate 
away from the wall. The overdamped character of the col- 
loidal dynamics is reflected by the absence of oscillations 
in the pressure as a function of time. 

The lifting of the COM after switch on changes the free 
energy, Eq. (fTQ)) , of the system, which, as the equilib- 
rium marks its minimum, has to increase. Its increase 
AF = AU - TAS = T(t) - T{t = 0), as shown in 
Fig. is nevertheless much smaller than the increase in 
potential energy V CKt , because the entropic energy —ST, 
given by the remaining terms in Eq. (fTX))) , decreases. AT 
shows most clearly the two timescales and the physical 
mechanisms at work: The initial increase up to roughly 
7T = 1 is due to the laning of particles near the wall, 
as observed in Fig. Q] Up to this point, the COM has 
hardly changed (compare with Fig. [2^,), and this initial 
increase is mostly due to decreasing entropy (see Fig. HJl). 



The entropy serves now as a driving force for the much 
slower migration of particles connected to the second in- 
crease in AJ-" to the final value, during which the potential 
energy increases and the entropic energy decreases, such 
that eventually AS(t — > 00) = 1.15 is positive. 

We note that after switch on of shear, energy is dis- 
sipated by two mechanisms: The vertical particle motion 
(in ^-direction), and the constant dissipation of energy due 
to the shear friction, connected to the shear stress, which 
persists even in the steady state. While the former can be 
derived using the particle currents (see Eq. (fT4")) below), 
the latter is at present not accessible in our theory, and 
we leave a full discussion of dissipation after switch on to 
future work. 

Given that the steady state has been reached, we now 
switch off the shear flow instantaneously at t = i ff (de- 
noting f = t — t Q ff) and investigate the relaxation back to 
equilibrium. Fig. [5t shows the subsequent decrease of the 
free energy back to the equilibrium value, which has to be 
dissipated via friction. The latter can be derived from the 
current of the one body density j(r, t) = p(y, t) ^ ^j^ffi 
in Eq. ([7]). Via Stokes' friction law, the force giving rise to 
this current is proportional to j, such that the following 
form arises for the energy dissipated after switch off 



df' j( - V ' T T = .F(f = 0)-.F(f). 



Wdiss(f) = / <Pr 



(14) 

The last equality can be shown analytically using the defi- 
nition in Eq. (|10|) and there is no need for numerical verifi- 
cation. Note that there is a difference between the current 



j and the average velocity 



dy(t) 
dt 



13]) : j depends on po- 



sition, and can be nonzero (due to relative motion) even 
if ^rP- is zero. 

dt 

We note that AS(f — > 00) < 0, as is typical for sed- 
imentation: The entropic energy — ST can be seen as a 
spring repelling the particles (or the COM) from the wall. 
As the COM sinks to lower values, the energy stored in 
the entropic spring naturally increases. 

In Fig[3] we focus on the time evolution of the density 
profile following switch off. For times shortly after ces- 
sation the heights of the peaks decrease (the main peak 
has its minimum value at f = 0.24 (see inset (a)), accom- 
panied with a reduction in the depth of the valleys, i.e., 
the peaks initially broaden, describing the inverse mech- 
anisms to the laning after switch on. This leads to the 
initial decrease in T seen in Fig. [2J because S initially in- 
creases. At long times the tail of the density distribution 
evolves monotonically, but now from the steady state back 
to equilibrium. 

A symmetry in the time evolution after switch on and 
off is apparent in our graphs. This can be understood by 
linearizing the density with respect to 7 after switch on 
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Fig. 3: Analogue to Fig[T] for shear cessation. The initial 
steady state (7 = 10) flow is switched off at f = 0, leading 
to a rapid transient 'undershoot' of the first peak followed by 
a slower approach to equilibrium. Inset (a) shows the height 
of the first peak as a function of time (in Brownian units). A 
minimum occurs for f = 0.24. Inset (b) shows the long range 
decay of the profile at f = and 00 (full black curves), and at 
intermediate times f = 2, 5 and 10 (broken red curves). 

and off, 

p(y,T)=p^(y)+jAp( yi T), (15) 
P(y, f) = p cq (y) + jAp(y, t ->• 00) - jAp(y, f ). (16) 

Of course, we have Ap{y, 00) = Ap{y, 00) such that 
the density approaches the equilibrium distribution after 
switch off. The straight forward linearization of Eq. ([7]) 
with respect to 7 yields the equality Ap(y,r) = Ap(y,f), 
stating the symmetry of the time evolution after switch 
on and off. From the appearance of this symmetry in our 
numerical results (compare the main panels if Figs. [T] and 
[3] as well as Fig.[2fc) we conclude that our chosen shear rate 
(despite 7 = 10 3> 1) is still not far away from this linear 
regime. We note that this symmetry does not depend on 
the details of the flow kernel in Eq. ([7]), and we expect it 
to be observable in experiments. However, as time grows, 
the full solution for finite shear rates deviates more and 
more from the linearized result, such that at long times, 
the density behaves in a more nontrivial way, leading to 
relaxation times which are not independent of shear rate 
(in contrast to the solutions for Ap(y,r)). 

Let us finally consider smooth variations in shear rate 
by applying oscillatory shear. By standard convention the 
strain in an oscillatory shear experiment is given by 

7(r) = 70 sin(wT)9(r), (17) 

where 70 is the strain amplitude, uj is the frequency and 
6(t) is the unit step function. The oscillatory shear pro- 
tocol (HH) can be realized experimentally, e.g. in a sim- 
ple parallel plate shear cell, which is particularly advanta- 
geous when seeking to combine rheological measurements 
with real space (confocal) microscopy (see e.g. [27]). 



In FiglU we show the time evolution of the centre-of- 
mass for four different frequencies to and strain amplitudes 
70, keeping the strain rate amplitude 70 u> (corresponding 
to the maximal shear rate) fixed at 70^ = 20. In all curves 
there is a transient dynamics of the COM, which initially 
follows the curve for steady shear with 7 = 20, consistent 
with linearization of the strain rate field for short times 
j(r) = 7owcos(o;t) w 7o^ ; until a periodic stationary 
state is attained. In this periodic state also the free energy 
(not shown) changes periodically while the system repeat- 
edly undergoes the mechanisms described above. For all 
curves shown, the stationary state is approached for times 
t f=s IOcj -1 , leading to an increasingly rapid deviation 
from the steady shear curve for increasing to. 

The stationary y(t) resembles a driven damped oscil- 
lator with frequency 2w (the response is invariant with 
respect to shear direction) and a phase shift 6 relative to 
the strain field (jTTJ) , such that its fundamental harmonic is 
~ sin(2o;r + 5). The inset of Fig. [4] shows S as a function 
of co, again for fixed 7oW- For uj — > 0, 7 changes slower 
and slower, and the COM y° sc in oscillatory flow can be 
expressed by its value y st (j) in steady shear, as 

lim y°-(r)=y st (7(r)). (18) 

u — >0,7o^J— const. 

In this limit, 5 thus approaches a value close to 7r/2, 
and the COM will oscillate between y st (j = 0) and 
y st {j = 7o^) (the equilibrium and steady shear lines in 
the figure). As to increases, 5 decreases (nearly exponen- 
tially), and the COM is more and more in phase with 
the strain. The phase shift S is reminiscent of the well 
known phase shift of the stress response with respect to 
strain [2j[TT], as described by the storage (G') and loss 
(G") moduli. G' measures the stress in phase with dis- 
tortion 7, giving the elastic response, while G" measures 
the stress in phase with rate 7, thereby quantifying viscous 
losses. As mentioned above, a detailed study of dissipation 
in our system (including COM motion) will be addressed 
in future 1 work, which appears indeed of importance for 
a complete description of viscous loss: If the colloids in 
an experiment are not perfectly density matched with the 
solvent (as is always the case in practice) the dissipation 
mechanism identified here, dragging the centre-of-mass up 
and down through the solvent, will provide an additional 
contribution to the loss modulus. 

As the frequency is increased, we observe that the os- 
cillations in y(r) decrease in amplitude. In this limit, the 
COM is too slow to follow the quick changes of the shear 
rate and the envelope function enclosing the oscillations 
converges smoothly to a finite value larger than the equi- 
librium value. It is intuitive that this value corresponds to 
the steady shear result taken at the average of I7I during 
one cycle, given by (I7I) = -^aH, We thus anticipate that 
the COM can also in this limit be expressed by its value 
in steady shear, 

lim ^( r )=y 8t f 7 = ^V (19) 

uj— *-oo,7oCl>= const. \ 71 J 
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Fig. 4: The time evolution of the centre-of-mass for the oscil- 
latory shear flow in Eq. (|17p for frequencies uj — 0.05 (brown), 
0.1 (blue), 0.5 (green) and 2 (red), with fixed 70W = 20. For 
long times a stationary state is attained. At high frequencies 
we sample only the average shear rate and the steady shear 
limit with 7 = 40/-7T = 12.73 is approached. The inset shows 
the phase shift of the COM relative to the applied shear strain. 
The calculated data points (circles) are consistent with an ex- 
ponential decay. 



In Fig. 2] we show the COM as function of time after 
switch on for 7 = 12.73, corresponding to the average 
rate for j q oj = 20. The present theory obeys for any 
fixed value of 70^. Physically, however, we would expect 
the COM to approach the equilibrium value for 70 -C 1, 
as the particles just move affinely back and forth over a 
strain range so small that they do not collide. We there- 
fore suspect that the flow kernel (J9j) only gives physically 
useful results for strains large enough that many particle 
collisions happen during any cycle. The fact that our ap- 
proach can only resolve time intervals within which many 
collisions have occurred represents an intrinsic limitation 
of our mean-field approach. 

Regarding an experimental sample, deviations from a 
linear shear profile, giving rise to deviations from Eq. (1191) . 
could be anticipated when lj^ 1 becomes comparable to 
the time taken for shear waves in the solvent to diffuse 
across the physical sample from one bounding plate to the 
other. Despite these limitations, the behavior described by 
Eq. (|19p should be observable in experiments for a certain 
frequency range. 

To conclude: Under time dependent shear the response 
of the density distribution of a colloidal sediment is also 
time dependent, as can be observed from its COM. After 
switch on of shear, a two step mechanism can be identi- 
fied, consisting of laning followed by migration. Several of 
the quantities discussed in the present work are accessible 
using real time particle tracking |27j and can be directly 
tested in experiment. Future work will focus on the time 
dependent shear stress and energy dissipation, aiming at 
a more complete description of the phenomenon of viscous 
resuspension. 
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